##R CODES
##
##
### BETA DISTRIBUTION
# 4P-Beta Distribution/LOGLIKELIHOOD FUNCTION		
lbeta4 <- function(theta,x,logyn=T,addsign=F){
# theta[1]=gamma_1, theta[2]=gamma_2, theta[3]=alpha_1 and theta[4]=alpha_2
n <- length(x)
res <- (-n*lbeta((theta[3]*(theta[3]>0)),(theta[4] *(theta[4]>0))))+(theta[3]-1)*sum(log((x-theta[1])*(x>=theta[1])))+(theta[4]-1)*sum(log((theta[2]-x)*(theta[2]>=x)))-n*(theta[3]+theta[4]-1)*log((theta[2]-theta[1])*(theta[2]>=theta[1]))
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}

# 2P-Beta Distribution/LOGLIKELIHOOD FUNCTION	
lbeta2 <- function(theta,x,logyn=T,addsign=F){
# theta[1]=alpha_1 and theta[2]=alpha_2
n <-length(x)
res <- (-n*lbeta(theta[1]*(theta[1]>0),theta[2]*(theta[2]>0)))+(theta[1]-1)*sum(log((x)*(x>=0)*(x<=1)))+(theta[2]-1)*sum(log((1-x)*(x>=0)*(x<=1)))
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}

# Beta Distribution/OPTIMIZATION FUNCTION
MLEbeta <- function(x,npar){
  # MLEs for 4P-Beta distribution:
  if (npar==4)  {

  # STARTING VALUES for 4P-Beta distribution:
  # Sample size:
  n <- length(x)

  # Boundary parameters; gamma_1 and gamma_2:	
  if (min(x)<0)  gamma_1 <- min(x)*(n+1)/(n)
  else           gamma_1 <- min(x)*(n)/(n+1)
  if (max(x)<0)  gamma_2 <- max(x)*(n)/(n+1)
  else           gamma_2 <- max(x)*(n+1)/(n)

  # Mean and variance are computed for 4P-beta distribution:
  a <- (mean(x)-gamma_1)/(gamma_2-gamma_1) # Mean
  v <- var(x)/(gamma_2-gamma_1)^2          # Variance

  # Shape parameters; alpha_1 and alpha_2:
  alpha_1 <- a*((a*(1-a)/v)-1)
  alpha_2 <- (1-a)*((a*(1-a)/v)-1)

  # We gather all starting values in the same vector:
  starting.beta <- c(gamma_1,gamma_2,alpha_1,alpha_2)

  # NUMERICAL SOLUTION for 4P-Beta distribution:
  res <- optim(starting.beta,lbeta4,x=x,addsign=T)
  r <- c(-2*(-res$value)+2*4,res$par)
  names(r) <- c("AIC","lower bound","upper bound","shape1","shape2")
  return (r)
  }
  
  # MLEs for 2P-Beta distribution:
  if (npar==2) {
  # distribution support for 2P-Beta distribution:
  if (min(x)<0|max(x)>1) return (NA)

  # STARTING VALUES for 2P-Beta distribution:

  # Mean and variance for 2P-beta distribution:
  a <- mean(x)      # Mean
  v <- var(x)       # Variance

  # Shape parameters; alpha_1 and alpha_2:
  alpha_1 <- a*((a*(1-a)/v)-1)
  alpha_2 <- (1-a)*((a*(1-a)/v)-1)
  
  # We gather all starting values in the same vector:
  starting.beta <- c(alpha_1,alpha_2)

  # NUMERICAL SOLUTION for 2P-Beta distribution
  res <- optim(starting.beta,lbeta2,x=x,addsign=T)
  r <- c(-2*(-res$value)+2*2,res$par)
  names(r) <- c("AIC","shape1","shape2")
  return (r)
  }
}

##### UNIFORM DISTRIBUTION
# 2P-Uniform Distribution/LOGLIKELIHOOD FUNCTION
luniform2 <- function(theta,x,logyn=T,addsign=F){
# theta[1]=a, theta[2]=b
n <- length(x)
res <- -n*(log((theta[2]-theta[1])*(theta[2]>=theta[1])))+sum(log(x>=theta[1]))+sum(log(x<=theta[2]))
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}

# 1P-Uniform Distribution/LOGLIKELIHOOD FUNCTION
luniform1 <- function(theta,x,logyn=T,addsign=F){
# theta=b
n <- length(x)
res <- -n*(log(theta))+sum(log(x>=0))+sum(log(x<=theta))
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}

# Uniform Distribution/OPTIMIZATION FUNCTION
MLEuniform <- function(x,npar){
  # MLEs for 2P-Uniform distribution:
  if (npar==2){
  # ANALYTICAL SOLUTION for 2P-Uniform distribution:
  # MLEs of boundary parameters; a and b:
  a_hat <- min(x)
  b_hat <- max(x)

  # LOG-LIKELHOOD VALUE FOR THE MLEs:
  res <- luniform2(c(a_hat,b_hat),x)
  r <- c(-2*(res)+2*2,a_hat,b_hat)
  names(r) <- c("AIC","lower bound","upper bound")
  return (r)
  }
  # MLE for 1P-Uniform distribution:
  if (npar==1){

  # distribution support for 2P-Uniform distribution:
  if (min(x)<0)  return (NA)

  # ANALYTICAL SOLUTION for 2P-Uniform distribution:
  # MLE of boundary parameter, b:
  b_hat <- max(x)

  # LOG-LIKELHOOD VALUE FOR THE MLE:
  res <- luniform1(b_hat,x)
  r <- c(-2*(res)+2*1,b_hat)
  names(r)<-c("AIC","upper bound")
  return (r)
  }
}

##### TRIANGULAR DISTRIBUTION
# Triangular Distribution/LOGLIKELIHOOD FUNCTION
ltriangular <- function(theta,x,logyn=T,addsign=F){
n <- length(x)
res <- sum((((log(2*(x-theta[1])*(x>=theta[1])))-log((theta[2]-theta[1])*(theta[2]>theta[1]))-log((theta[3]-theta[1])*(theta[3]>=theta[1])))*((x>=theta[1])*(x<=theta[3])==1))+(((log(2*(theta[2]-x)*(x<=theta[2])))-log((theta[2]-theta[1])*(theta[2]>theta[1]))-log((theta[2]-theta[3])*(theta[3]<=theta[2])))*((x>theta[3])*(x<=theta[2])==1)))
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}
# Triangular Distribution/OPTIMIZATION FUNCTION
MLEtriangular <- function(x){

  # STARTING VALUES for Triangular distribution:
  # Sample size:
  n <- length(x)

  # Boundary parameters; a and b:
  if  (min(x)<0)  a <- min(x)*(n+1)/(n)
  else            a <- min(x)*(n)/(n+1)
  if  (max(x)<0)  b <- max(x)*(n)/(n+1)
  else            b <- max(x)*(n+1)/n
 
  # Shape parameter (mode value); c:
  Mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
  }
  c <- Mode(x) 

  # We gather all starting values in the same vector:
  starting.triangular <- c(a,b,c)

  # NUMERICAL SOLUTION for Triangular distribution:
  res <- optim(starting.triangular,ltriangular,x=x,addsign=T)
  r <- c(-2*(-res$value)+2*3,res$par)
  names(r) <- c("AIC","lower bound","upper bound","mode")
  return (r)
}

##### GAMMA DISTRIBUTION
# 3P-Gamma Distribution/PROFILE LOGLIKELIHOOD FUNCTION
lgamma3 <- function(theta,x,logyn=T,addsign=F){
#theta[1]=alpha, theta[2]=gamma
n <- length(x)
res <- -n*(theta[1]*(theta[1]>0))*log(sum((x-theta[2])*(x>=theta[2])))+n*(theta[1]*(theta[1]>0))*log((n*theta[1]*(theta[1]>0)))+(theta[1]*(theta[1]>0)-1)*sum(log((x-theta[2])*(x>=theta[2])))-(n*theta[1]*(theta[1]>0))-n*lgamma(theta[1]*(theta[1]>0))
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}
# 2P-Gamma Distribution/PROFILE LOGLIKELIHOOD FUNCTION
lgamma2 <- function(theta,x,logyn=T,addsign=F){
#theta[1]=alpha
n <- length(x)
res <- -n*(theta*(theta>0))*log(sum((x)*(x>=0)))+n*(theta*(theta>0))*log((n*theta*(theta>0)))+(theta*(theta>0)-1)*sum(log((x)*(x>=0)))-(n*theta*(theta>0))-n*lgamma(theta*(theta>0))
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}
# Gamma Distribution/OPTIMIZATION FUNCTION
MLEgamma<-function(x,npar){
  # MLEs for 3P-Gamma distribution:
  if (npar==3) {

  # STARTING VALUES for 3P-Gamma distribution:
  # Sample size:
  n <- length(x)

  # Location parameter, gamma:
  if (min(x)<0)  gamma <- min(x)*(n+1)/(n)   
  else           gamma <- min(x)*(n)/(n+1)
  
  # Shape parameter, alpha:
  xs <- x-gamma
  alpha <- (mean(xs)/sd(xs))^2

  # We gather all starting values in the same vector:
  starting.gamma <- c(alpha,gamma)

  # NUMERICAL SOLUTION for 3P-Gamma distribution:
  res <- optim(starting.gamma,lgamma3,x=x,addsign=T)
  
  # Scale parameter, beta:
  beta_hat <- sum(x-res$par[2])/(n*res$par[1])

  r <- c(-2*(-res$value)+2*3,c(res$par,beta_hat))
  names(r)<-c("AIC","shape","location","scale")
  return (r)
  }

  # MLEs for 2P-Gamma distribution:
  if (npar==2){
  # distribution support for 2P-Gamma distribution:
  if(min(x)<0)     return(NA)

  # STARTING VALUE for 2P-Gamma distribution:
  # Shape parameter, alpha:
  alpha <- (mean(x)/sd(x))^2
  starting.gamma <- alpha

  # NUMERICAL SOLUTION for 2P-Gamma distribution:
  res <- optim(starting.gamma,lgamma2,method="BFGS",x=x,addsign=T)
  
  # Scale parameter, beta:
  n <- length(x)
  beta <- sum(x)/(n*res$par)

  r <- c(-2*(-res$value)+2*2,c(res$par,beta))
  names(r) <- c("AIC","shape","scale")
  return (r)
  }
}

##### ERLANG DISTRIBUTION
# Erlang Distribution/PROFILE LOGLIKELIHOOD FUNCTION
lerlang2 <- function(theta,x,logyn=T,addsign=F){
# theta=alpha
n <- length(x)
res <- -n*as.integer(theta*(theta>0))*log(sum((x)*(x>=0)))+n*as.integer(theta*(theta>0))*log((n*as.integer(theta*(theta>0))))+as.integer(theta*(theta>0)-1)*sum(log((x)*(x>=0)))-(n*as.integer(theta*(theta>0)))-n*lgamma(as.integer(theta*(theta>0)))
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}

# Erlang Distribution/OPTIMIZATION FUNCTION
MLEerlang<-function(x){
  # MLEs for Erlang distribution:

  # distribution support for 2P-Erlang distribution:
  if(min(x)<0)     return(NA)


  # STARTING VALUE for 2P-Erlang distribution:
  # Shape parameter, alpha:
  k <- MLEgamma(x,npar=3)[2]
  if (min(x)<0) {alpha <- NA}
  else if (k<1) {alpha <- 1}
  else          {alpha <- round(k)}
  starting.erlang <- alpha
    
  # NUMERICAL SOLUTION for 2P-Erlang distribution:
  res <- optim(starting.erlang,lerlang2,method="BFGS",x=x,addsign=T)
  
  # Scale parameter, beta:
  n <- length(x)
  beta <- sum(x)/(n*res$par)

  r <- c(-2*(-res$value)+2*2,c(res$par,beta))
  names(r) <- c("AIC","shape","scale")
  return (r)
 }

##### WEIBULL DISTRIBUTION
# 3P-Weibull Distribution/PROFILE LOGLIKELIHOOD FUNCTION
lweibull3 <- function(theta,x,logyn=T,addsign=F){
# theta[1]=alpha, theta[2]=gamma
n <- length(x)
res <- n*log(theta[1]*(theta[1]>0))-n*theta[1]*log(((sum((x-theta[2])^(theta[1])))/n)^(1/theta[1]))+(theta[1]-1)*sum(log((x-theta[2])*(x>=theta[2])))-sum(((x-theta[2])*(x>=theta[2])/((sum((x-theta[2])^(theta[1])))/n)^(1/theta[1]))^theta[1])
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}

# 2P-Weibull Distribution/PROFILE LOGLIKELIHOOD FUNCTION
lweibull2 <- function(theta,x,logyn=T,addsign=F){
# theta=alpha
n <- length(x)
res <- n*log(theta*(theta>0))-n*theta*log(((sum(x^(theta)))/n)^(1/theta))+(theta-1)*sum(log((x)*(x>=0)))-sum(((x)*(x>=0)/((sum(x^(theta)))/n)^(1/theta))^theta)
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}

# Weibull Distribution/OPTIMIZATION FUNCTION
MLEweibull <- function(x,npar){
  # MLEs for 3P-Weibull distribution:
  if (npar==3){
  # STARTING VALUES for 3P-Weibull distribution:
  # Location parameter, gamma:
  if (min(x)<0)  gamma <- min(x)+0.1*min(x)  
  else           gamma <- min(x)-0.1*min(x)

  # Shape parameter, alpha:
  xs <- x-gamma
  alpha <- mean(xs)/sd(xs)

  # We gather all starting values in the same vector:
  starting.weibull <- c(alpha,gamma)

  # NUMERICAL SOLUTION for 3P-Weibull distribution:
  res <- optim(starting.weibull,lweibull3,x=x,addsign=T)
  n <- length(x)

  # Scale parameter, beta:
  beta_hat <- ((sum((x-res$par[2])^(res$par[1])))/n)^(1/res$par[1])
 
  r <- c(-2*(-res$value)+2*3,c(res$par,beta_hat))
  names(r) <- c("AIC","shape","location","scale")
  return (r)
  }
  # MLEs for 2P-Weibull distribution:
  if (npar==2){
  # distribution support for 2P-Weibull distribution:
  if(min(x)<0)     return(NA)

  # STARTING VALUE for 2P-Weibull distribution:
  # Shape parameter=alpha:
  alpha <- mean(x)/sd(x)
  if (alpha<0) return (NA)
  starting.weibull <- alpha

  # NUMERICAL SOLUTION for 2P-Weibull distribution:
  # if(!is.finite(lweibull2(starting.weibull,x))) return(res<-NA)
  res <- optim(starting.weibull,lweibull2,method="BFGS",x=x,addsign=T)

  # Scale parameter, beta:
  n <- length(x)
  beta_hat <- ((sum(x^(res$par)))/n)^(1/res$par)

  r <- c(-2*(-res$value)+2*2,c(res$par,beta_hat))
  names(r) <- c("AIC","shape","scale")
  return (r)
  }
}

##### EXPONENTIAL DISTRIBUTION
# 2P-Exponential Distribution/LOGLIKELIHOOD FUNCTION
lexponential2 <- function(theta,x,logyn=T,addsign=F){
# theta[1]=beta, theta[2]=gamma
n <- length(x)
res <- n*log(theta[1]*(theta[1]>0))-((theta[1])*(theta[1]>0)*sum((x-theta[2])*(x>=theta[2])))+sum(log(x>=theta[2]))
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}
# 1P-Exponential Distribution/LOGLIKELIHOOD FUNCTION
lexponential1 <- function(theta,x,logyn=T,addsign=F){
# theta=beta
n <- length(x)
res <- n*log(theta*(theta>0))-((theta)*(theta>0)*sum((x)*(x>=0)))+sum(log(x>=0))
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}

# Exponential Distribution/OPTIMIZATION FUNCTION
MLEexponential <- function(x,npar){
  # MLEs for 2P-Exponential distribution:
  if(npar==2){
  # ANALYTICAL SOLUTION for 2P-Exponential distribution:
  # MLEs of location and rate parameters:
    # location parameter, gamma:
    gamma_hat <- min(x)

    # rate parameter, beta:
    beta_hat <- 1/mean(x-gamma_hat)

  # LOG-LIKELHOOD VALUE FOR THE MLEs:
   res <- lexponential2(c(beta_hat,gamma_hat),x)
   r <- c(-2*(res)+2*2,beta_hat,gamma_hat)
   names(r) <- c("AIC","rate","location")
   return (r)
   }
  # MLEs for 1P-Exponential distribution:
   if(npar==1){
  # distribution support for 1P-Exponential distribution:
   if (min(x)<0)  return (NA)

  # MLEs of rate parameter, beta:
   beta_hat <- 1/mean(x)  

  # LOG-LIKELHOOD VALUE FOR THE MLE:
   res <- lexponential1(beta_hat,x)
   r <- c(-2*(res)+2*1,beta_hat)
   names(r)<-c("AIC","rate")
   return (r)
   }
}

##### NORMAL DISTRIBUTION
# Normal Distribution/LOGLIKELIHOOD FUNCTION
lnormal <- function(theta,x,logyn=T,addsign=F){
# theta[1]=mu,theta[2]=sigma
n <- length(x)
res <- -n*log(sqrt(2*pi))-n*log(theta[2]*(theta[2]>0))-sum((((x)-theta[1])^2)/(2*(theta[2]*(theta[2]>0))^2))
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}

# Normal Distribution/OPTIMIZATION FUNCTION
MLEnormal <- function(x){
  # MLEs for Normal distribution:
  # ANALYTICAL SOLUTION for Normal distribution:
  # MLEs of location and scale parameters:
    # Location parameter, mu:
    n <- length(x)
    mu_hat <- mean(x)

    # Scale parameter, sigma: 
    sigma_hat <- sqrt((1/n)*sum((x-mu_hat)^2))

  # LOG-LIKELHOOD VALUE FOR THE MLEs:    
    res <- lnormal(c(mu_hat,sigma_hat),x)
    r <- c(-2*(res)+2*2,mu_hat,sigma_hat)
    names(r) <- c("AIC","location","scale")
    return (r)
}

##### STUDENT'S T-DISTRIBUTION
# 3P-Student's T Distribution/LOGLIKELIHOOD FUNCTION	
lt3 <- function(theta,x,logyn=T,addsign=F){
# theta[1]=mu,theta[2]=sigma,theta[3]=nu(degrees of freedom)
n <- length(x)
res <- n*lgamma((theta[3]*(theta[3]>0)+1)/2)-n*lgamma(theta[3]*(theta[3]>0)/2)-n*log(theta[2]*(theta[2]>0))-n*(1/2)*log(pi*theta[3]*(theta[3]>0))-((theta[3]*(theta[3]>0)+1)/2)*sum(log(1+(((x-theta[1])/theta[2])^2)/(theta[3]*(theta[3]>0))))
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}

# 1P-Student's T Distribution/LOGLIKELIHOOD FUNCTION	
lt1 <- function(theta,x,logyn=T,addsign=F){
# theta=nu(degrees of freedom)
n <- length(x)
res <- (-n/2)*log(theta*(theta>0))-n*lbeta((1/2),((theta*(theta>0))/2))-((theta*(theta>0)+1)/2)*sum(log(1+((x)^2)/theta*(theta>0)))
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}

MLEt <- function(x,npar){
  # MLEs for 3P-Student's T distribution:
  if (npar==3){

  # STARTING VALUES for 3P-Student's T Distribution:
  # Location parameter, gamma:
  gamma <- mean(x)
  
  # Scale parameter, sigma:
  sigma <- 0.01

  # Degrees of freedom, nu:
  nu <- 0.01

  # We gather all starting values in the same vector:
  starting.st <- c(gamma,sigma,nu) 

  # NUMERICAL SOLUTION for 3P-Student's T Distribution:
  res <- optim(starting.st,lt3,x=x,addsign=T)
  r <- c(-2*(-res$value)+2*3, res$par)
  names(r) <- c("AIC","location","scale","degrees of freedom")
  return (r)
  }

  # MLEs for 1P-Student's T distribution:
  if (npar==1){

  # STARTING VALUE for 1P-Student's T Distribution:
  # Degrees of freedom, nu:
  nu_lower <- 0.01
  nu_upper <- 500

  starting.st <- c(nu_lower,nu_upper) 
  
  # NUMERICAL SOLUTION for 1P-Student's T Distribution:
  res <- optimize(lt1, x=x, starting.st, maximum=TRUE)
  r <- c(-2*(res$objective)+2*1, res$maximum)
  names(r) <- c("AIC","degrees of freedom")
  return (r)
  }
}

###### LOGNORMAL DISTRIBUTION
# 3P-Lognormal Distribution/PROFILE LOGLIKELIHOOD FUNCTION
llognormal3 <- function(theta,x,logyn=T,addsign=F){
# theta=gamma
n <- length(x)
res <- -(n*(log(sqrt(2*pi))))-(n*log(sqrt(sum((log((x-theta)*(x>theta))-sum(log((x-theta)*(x>theta))/n))^2)/n)))-sum(log((x-theta)*(x>theta)))-sum((((log(((x-theta))*(x>theta)))-sum(log((x-theta)*(x>theta))/n))^2)/(2*(sqrt(sum((log((x-theta)*(x>theta))-sum(log((x-theta)*(x>theta))/n))^2)/n))^2))
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}
# 2P-Lognormal Distribution/LOGLIKELIHOOD FUNCTION
llognormal2 <- function(theta,x,logyn=T,addsign=F){
# theta[1]=mu and theta[2]=sigma
n <- length(x)
res <- -(n*(log(sqrt(2*pi))))-(n*log(theta[2]*(theta[2]>0)))-sum(log(x*(x>0)))-sum((((log((x)*(x>0)))-theta[1])^2)/(2*(theta[2]*(theta[2]>0))^2))
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}

# Lognormal Distribution/OPTIMIZATION FUNCTION
MLElognormal <- function(x,npar){
  # MLEs for 3P-Lognormal distribution:
  if(npar==3){

  # STARTING VALUES for 3P-Lognormal distribution:
  n <- length(x)
  # Location parameter, gamma:
  if (min(x)<0)   gamma <- min(x)*(n+1)/(n)
  else            gamma <- min(x)*(n)/(n+1)
 
  gamma_upper <- gamma
  gamma_lower <- -2000000 #lower limit for gamma which is a dummy number

  starting.lognormal <- c(gamma_lower,gamma_upper) 

  # NUMERICAL SOLUTION for 3P-Lognormal distribution:
  res <- optimize(llognormal3,x=x,starting.lognormal,addsign=T)
  
  gamma_hat <- res$minimum

  # Scale parameter, mu:
  mu_hat <- sum(log((x-gamma_hat))/n)

  # Shape parameter, sigma:
  sigma_hat <- sqrt(sum((log(x-gamma_hat)-mu_hat)^2)/n)

  r <- c(-2*(-res$objective)+2*3,mu_hat,sigma_hat,gamma_hat)
  names(r) <- c("AIC","scale","shape","location")
  return (r)
  }
  # MLEs for 2P-Lognormal distribution:
  if(npar==2){
  # distribution support for 2P-Lognormal distribution:
  if (min(x)<=0)    return (NA)
 
  # ANALYTICAL SOLUTION for 2P-Lognormal distribution
  # MLEs of scale and shape parameters:
    # Scale parameter, mu:
    n <- length(x)
    mu_hat <- sum(log(x)/n)

    # Shape parameter, sigma: 
    sigma_hat <- sqrt(sum((log(x)-mu_hat)^2)/n)

  # LOG-LIKELHOOD VALUE FOR THE MLEs: 
  res <- llognormal2(c(mu_hat,sigma_hat),x)
  r <- c(-2*(res)+2*2,mu_hat,sigma_hat)
  names(r) <- c("AIC","scale","shape")
  return (r)
  }
}



##### LOGISTIC DISTRIBUTION
# Logistic Distribution/LOGLIKELIHOOD FUNCTION
llogistic <- function(theta,x,logyn=T,addsign=F){
# theta[1]=mu and theta[2]=alpha
n <- length(x)
res <- (-sum(x-theta[1])/(theta[2]*(theta[2]>0)))-n*log((theta[2]*(theta[2]>0)))-2*sum(log(1+exp(-(x-theta[1])/(theta[2]*(theta[2]>0)))))
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}
# Logistic Distribution/OPTIMIZATION FUNCTION
MLElogistic <- function(x){

  # MLEs for Logistic distribution:
  # STARTING VALUES for Logistic distribution:
  n <- length(x)
  # Location parameter, mu:
  mu <- mean(x)
 
  # Scale parameter, alpha:
  alpha <- (sqrt(3)/pi)*((sum(x^2)/n)-mu^2)^(1/2)

  # We gather all starting values in the same vector:
  starting.logistic <- c(mu,alpha) 
 
  # NUMERICAL SOLUTION for Logistic distribution:
  # if(!is.finite(llogistic(starting.logistic,x))) return(NA)

  res <- optim(starting.logistic,llogistic,x=x,addsign=T)
  r <- c(-2*(-res$value)+2*2, res$par)
  names(r) <- c("AIC","location","scale")
  return (r)
}

#####RAYLEIGH DISTRIBUTION
# 2P-Rayleigh Distribution/PROFILE LOGLIKELIHOOD FUNCTION	
lrayleigh2 <- function(theta,x,logyn=T,addsign=F){
# theta=gamma
n <- length(x)
res <- sum(log((x-theta)*(x>=theta)))-2*n*log(sqrt((sum(((x-theta)*(x>=theta))^2))/(2*n)))-(sum(((x-theta)*(x>=theta))^2)/(2*(sqrt((sum(((x-theta)*(x>=theta))^2))/(2*n)))^2))
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}

# 1P-Rayleigh Distribution/LOGLIKELIHOOD FUNCTION
lrayleigh1 <- function(theta,x,logyn=T,addsign=F){
# theta=sigma
n <- length(x)
res <- sum(log(x*(x>=0)))-2*n*log((theta)*(theta>0))-(sum((x*(x>=0))^2)/(2*(theta*(theta>0))^2))
if (!logyn)  res <- exp(res)
if (addsign) res <- -res
res
}
# Rayleigh Distribution/OPTIMIZATION FUNCTION
MLErayleigh <- function(x,npar){
  # MLEs for 2P-Rayleigh distribution:
  if (npar==2){

  # STARTING VALUES for 2P-Rayleigh distribution:
  n <- length(x)
  # Location parameter, gamma:
  if (min(x)<0) gamma <- min(x)*(n+1)/(n) 
  else          gamma <- min(x)*(n)/(n+1)

  gamma_upper <- gamma
  gamma_lower <- -20000000 #lower limit for gamma which is a dummy number

  starting.rayleigh <- c(gamma_lower,gamma_upper) 

  # NUMERICAL SOLUTION for 2P-Rayleigh distribution:
  res <- optimize(lrayleigh2,x=x,starting.rayleigh,addsign=T)

  gamma_hat <- res$minimum
  # Scale parameter, sigma:
  sigma_hat <- sqrt((sum((x-gamma_hat)^2))/(2*n))

  r <- c(-2*(-res$objective)+2*2,gamma_hat,sigma_hat)
  names(r) <- c("AIC","location","scale")
  return (r)
  }
  # MLE for 1P-Rayleigh distribution:
  if (npar==1){
  # distribution support for 1P-Rayleigh distribution:
  if (min(x)<0) return (NA)

  # ANALYTICAL SOLUTION for 1P-Rayleigh distribution:
  n <- length(x)
  # Scale parameter, sigma:
  sigma_hat <- sqrt((sum(x^2))/(2*n))


  # LOG-LIKELHOOD VALUE FOR THE MLE: 
  res <- (lrayleigh1(sigma_hat,x))
  r <- c(-2*(res)+2*1,sigma_hat)
  names(r) <- c("AIC","scale")
  return (r)
  }
}


#library(VGAM)#vglm(x ~ 1, rayleigh, trace=TRUE)#Coef(fit)

alldistributions<-function(x){
result<-matrix(NA,nrow=25,ncol=5)
rownames(result)<-c("Beta(4)","Beta(2)","Uniform(2)","Uniform(1)","Triangular","Gamma(3)","Gamma(2)","Erlang(2)","Weibull(3)","Weibull(2)","Exponential(2)","Exponential(1)","Normal","Student's-t(3)","Student's-t(1)","Lognormal(3)","Lognormal(2)","Logistic","Rayieigh(2)","Rayleigh(1)","N-Gamma(3)","N-Weibull(3)","N-Exponential(2)","N-Lognormal(3)","N-Rayleigh(2)")
colnames(result)<-c("AIC","p1","p2","p3","p4")              
result[1,1:5]<-MLEbeta(x,npar=4)
result[2,1:3]<-MLEbeta(x,npar=2)
result[3,1:3]<-MLEuniform(x,npar=2)
result[4,1:2]<-MLEuniform(x,npar=1)
result[5,1:4]<-MLEtriangular(x)
result[6,1:4]<-MLEgamma(x,npar=3)
result[7,1:3]<-MLEgamma(x,npar=2)
result[8,1:3]<-MLEerlang(x)
result[9,1:4]<-MLEweibull(x,npar=3)
result[10,1:3]<-MLEweibull(x,npar=2)
result[11,1:3]<-MLEexponential(x,npar=2)
result[12,1:2]<-MLEexponential(x,npar=1)
result[13,1:3]<-MLEnormal(x)
result[14,1:4]<-MLEt(x,npar=3)
result[15,1:2]<-MLEt(x,npar=1)
result[16,1:4]<-MLElognormal(x,npar=3)
result[17,1:3]<-MLElognormal(x,npar=2)
result[18,1:3]<-MLElogistic(x)
result[19,1:3]<-MLErayleigh(x,npar=2)
result[20,1:2]<-MLErayleigh(x,npar=1)
result[21,1:4]<-MLEgamma(-x,npar=3)
result[22,1:4]<-MLEweibull(-x,npar=3)
result[23,1:3]<-MLEexponential(-x,npar=2)
result[24,1:4]<-MLElognormal(-x,npar=3)
result[25,1:3]<-MLErayleigh(-x,npar=2)
result
}

par<-function(result){
result<-matrix(NA,nrow=25,ncol=4)
rownames(result)<-c("Beta(4)","Beta(2)","Uniform(2)","Uniform(1)","Triangular","Gamma(3)","Gamma(2)","Erlang(2)","Weibull(3)","Weibull(2)","Exponential(2)","Exponential(1)","Normal","Student's-t(3)","Student's-t(1)","Lognormal(3)","Lognormal(2)","Logistic","Rayieigh(2)","Rayleigh(1)","N-Gamma(3)","N-Weibull(3)","N-Exponential(2)","N-Lognormal(3)","N-Rayleigh(2)")
colnames(result)<-c("pname1","pname2","pname3","pname4")              
result[1,1:4]<-c("lower bound","upper bound","shape1","shape2")
result[2,1:2]<-c("shape1","shape2")
result[3,1:2]<-c("lower bound","upper bound")
result[4,1:1]<-"upper bound"
result[5,1:3]<-c("lower bound","upper bound","mode")
result[6,1:3]<-c("shape","location","scale")
result[7,1:2]<-c("shape","scale") 
result[8,1:2]<-c("shape","scale") 
result[9,1:3]<-c("shape","location","scale")
result[10,1:2]<-c("shape","scale") 
result[11,1:2]<-c("rate","location")
result[12,1:1]<-"rate"
result[13,1:2]<-c("location=mean","scale=sdev")
result[14,1:3]<-c("location","scale","degrees of freedom")
result[15,1:1]<-"degrees of freedom"
result[16,1:3]<-c("meanlog","sdlog","location")
result[17,1:2]<-c("meanlog","sdlog")
result[18,1:2]<-c("location","scale")
result[19,1:2]<-c("location","scale")
result[20,1:1]<-"scale"
result[21,1:3]<-c("shape","location","scale")
result[22,1:3]<-c("shape","location","scale")
result[23,1:2]<-c("rate","location")
result[24,1:3]<-c("meanlog","sdlog","location")
result[25,1:2]<-c("location","scale")
result
}

parnames<-par(result)

 m <-alldistributions(x)
 ms<-m[order(m[,1]),]
 ms[!is.na(ms[,1]),]

# CHI-SQUARE TEST FOR ALL DISTRIBUTIONS:
cs.test<-function(x){
options(warn=-1)
# PARAMETER ESTIMATION FOR ALL DISTRIBUTIONS:
  m <- alldistributions(x)

# SAMPLE SIZE:
  n <- length(x)  

# PREVIOUSLY UNDEFINED CDFS IN R: 
ptriangle1<-function(x,a,b,c){
res<- (((x-a)^2)/((b-a)*(c-a)))*(x>=a)*(x<=c)+(1-((b-x)^2/((b-a)*(b-c))))*(b>=x)*(x>c)+0*(x<a)+0*(x>b)
res
}
prayleigh1<-function(x,a){
res<-1-exp(-(x*(x>=0))^2/(2*(a*(a>0))^2))
res
}

# NUMBER OF CLASS INTERVALS:
k<-as.integer(sqrt(length(x)))

# CONVERT UNIFORM DENSITY:          
unif.beta4        <- pbeta((x-m[1,2])/(m[1,3]-m[1,2]),m[1,4],m[1,5])
unif.beta2        <- pbeta(x,m[2,2],m[2,3])
unif.uniform2     <- punif(x,m[3,2],m[3,3])
unif.uniform1     <- punif(x,0,m[4,2])
unif.triangular   <- ptriangle1(x,m[5,2],m[5,3],m[5,4])
unif.gamma3       <- pgamma(x-m[6,3],m[6,2],scale=m[6,4])
unif.gamma2       <- pgamma(x,m[7,2],scale=m[7,3])
unif.erlang2      <- pgamma(x,m[8,2],scale=m[8,3])
unif.weibull3     <- pweibull(x-m[9,3],m[9,2],scale=m[9,4])
unif.weibull2     <- pweibull(x,m[10,2],scale=m[10,3])
unif.exponential2 <- pexp(x-m[11,3],rate=m[11,2])
unif.exponential1 <- pexp(x,rate=m[12,2])
unif.normal       <- pnorm(x,mean=m[13,2],sd=m[13,3])
unif.t3           <- pt(((x-m[14,2])/m[14,3]),m[14,4])
unif.t1           <- pt(x,m[15,2])
unif.lognormal3   <- plnorm(x-m[16,4],m[16,2],m[16,3])
unif.lognormal2   <- plnorm(x,m[17,2],m[17,3])
unif.logistic     <- plogis(x,m[18,2],m[18,3])
unif.rayleigh2    <- prayleigh1(x-m[19,2],m[19,3])
unif.rayleigh1    <- prayleigh1(x,m[20,2])
unif.ngamma3      <- pgamma(-x-m[21,3],m[21,2],scale=m[21,4])
unif.nweibull3    <- pweibull(-x-m[22,3],m[22,2],scale=m[22,4])
unif.nexponential2<- pexp(-x-m[23,3],rate=m[23,2])
unif.nlognormal3  <- plnorm(-x-m[24,4],m[24,2],m[24,3])
unif.nrayleigh2   <- prayleigh1(-x-m[25,2],m[25,3])

# FIND OBSERVED FREQUENCY FOR EACH CLASS INTERVAL:
obs.beta4         <- table(floor(unif.beta4*k))
obs.beta2         <- table(floor(unif.beta2*k))   
obs.uniform2      <- table(floor(unif.uniform2*k))
obs.uniform1      <- table(floor(unif.uniform1*k))
obs.triangular    <- table(floor(unif.triangular*k))
obs.gamma3        <- table(floor(unif.gamma3*k))
obs.gamma2        <- table(floor(unif.gamma2*k))
obs.erlang2       <- table(floor(unif.erlang2*k))
obs.weibull3      <- table(floor(unif.weibull3*k))
obs.weibull2      <- table(floor(unif.weibull2*k))
obs.exponential2  <- table(floor(unif.exponential2*k))
obs.exponential1  <- table(floor(unif.exponential1*k))
obs.normal        <- table(floor(unif.normal*k))
obs.t3            <- table(floor(unif.t3*k))
obs.t1            <- table(floor(unif.t1*k))
obs.lognormal3    <- table(floor(unif.lognormal3*k))
obs.lognormal2    <- table(floor(unif.lognormal2*k))
obs.logistic      <- table(floor(unif.logistic*k))
obs.rayleigh2     <- table(floor(unif.rayleigh2*k))
obs.rayleigh1     <- table(floor(unif.rayleigh1*k))
obs.ngamma3       <- table(floor(unif.ngamma3*k))
obs.nweibull3     <- table(floor(unif.nweibull3*k))
obs.nexponential2 <- table(floor(unif.nexponential2*k))
obs.nlognormal3   <- table(floor(unif.nlognormal3*k))
obs.nrayleigh2    <- table(floor(unif.nrayleigh2*k))

# CHI-SQUARE TEST STATISTIC FOR EACH DISTRIBUTION:
if(length(obs.beta4)         >2 & is.finite(sum(unif.beta4)))         ct.beta4         <-(chisq.test(obs.beta4,        p = rep(1/length(obs.beta4),        length(obs.beta4))))        $statistic  else  ct.beta4        <-NA
if(length(obs.beta2)         >2 & is.finite(sum(unif.beta2)))         ct.beta2         <-(chisq.test(obs.beta2,        p = rep(1/length(obs.beta2),        length(obs.beta2))))        $statistic  else  ct.beta2        <-NA
if(length(obs.uniform2)      >2 & is.finite(sum(unif.uniform2)))      ct.uniform2      <-(chisq.test(obs.uniform2,     p = rep(1/length(obs.uniform2),     length(obs.uniform2))))     $statistic  else  ct.uniform2     <-NA
if(length(obs.uniform1)      >2 & is.finite(sum(unif.uniform1)))      ct.uniform1      <-(chisq.test(obs.uniform1 ,    p = rep(1/length(obs.uniform1),     length(obs.uniform1))))     $statistic  else  ct.uniform1     <-NA
if(length(obs.triangular)    >2 & is.finite(sum(unif.triangular)))    ct.triangular    <-(chisq.test(obs.triangular,   p = rep(1/length(obs.triangular),   length(obs.triangular))))   $statistic  else  ct.triangular   <-NA
if(length(obs.gamma3)        >2 & is.finite(sum(unif.gamma3)))        ct.gamma3        <-(chisq.test(obs.gamma3,       p = rep(1/length(obs.gamma3),       length(obs.gamma3))))       $statistic  else  ct.gamma3       <-NA
if(length(obs.gamma2)        >2 & is.finite(sum(unif.gamma2)))        ct.gamma2        <-(chisq.test(obs.gamma2,       p = rep(1/length(obs.gamma2),       length(obs.gamma2))))       $statistic  else  ct.gamma2       <-NA
if(length(obs.erlang2)       >2 & is.finite(sum(unif.erlang2)))       ct.erlang2       <-(chisq.test(obs.erlang2,      p = rep(1/length(obs.erlang2),      length(obs.erlang2))))      $statistic  else  ct.erlang2      <-NA
if(length(obs.weibull3)      >2 & is.finite(sum(unif.weibull3)))      ct.weibull3      <-(chisq.test(obs.weibull3,     p = rep(1/length(obs.weibull3),     length(obs.weibull3))))     $statistic  else  ct.weibull3     <-NA
if(length(obs.weibull2)      >2 & is.finite(sum(unif.weibull2)))      ct.weibull2      <-(chisq.test(obs.weibull2,     p = rep(1/length(obs.weibull2),     length(obs.weibull2))))     $statistic  else  ct.weibull2     <-NA
if(length(obs.exponential2)  >2 & is.finite(sum(unif.exponential2)))  ct.exponential2  <-(chisq.test(obs.exponential2, p = rep(1/length(obs.exponential2), length(obs.exponential2)))) $statistic  else  ct.exponential2 <-NA
if(length(obs.exponential1)  >2 & is.finite(sum(unif.exponential1)))  ct.exponential1  <-(chisq.test(obs.exponential1, p = rep(1/length(obs.exponential1), length(obs.exponential1)))) $statistic  else  ct.exponential1 <-NA
if(length(obs.normal)        >2 & is.finite(sum(unif.normal)))        ct.normal        <-(chisq.test(obs.normal,       p = rep(1/length(obs.normal),       length(obs.normal))))       $statistic  else  ct.normal       <-NA
if(length(obs.t3)            >2 & is.finite(sum(unif.t3)))            ct.t3            <-(chisq.test(obs.t3,           p = rep(1/length(obs.t3),           length(obs.t3))))           $statistic  else  ct.t3           <-NA
if(length(obs.t1)            >2 & is.finite(sum(unif.t1)))            ct.t1            <-(chisq.test(obs.t1,           p = rep(1/length(obs.t1),           length(obs.t1))))           $statistic  else  ct.t1           <-NA
if(length(obs.lognormal3)    >2 & is.finite(sum(unif.lognormal3)))    ct.lognormal3    <-(chisq.test(obs.lognormal3,   p = rep(1/length(obs.lognormal3),   length(obs.lognormal3))))   $statistic  else  ct.lognormal3   <-NA
if(length(obs.lognormal2)    >2 & is.finite(sum(unif.lognormal2)))    ct.lognormal2    <-(chisq.test(obs.lognormal2,   p = rep(1/length(obs.lognormal2),   length(obs.lognormal2))))   $statistic  else  ct.lognormal2   <-NA
if(length(obs.logistic)      >2 & is.finite(sum(unif.gamma3)))        ct.logistic      <-(chisq.test(obs.logistic,     p = rep(1/length(obs.logistic),     length(obs.logistic))))     $statistic  else  ct.logistic     <-NA
if(length(obs.rayleigh2)     >2 & is.finite(sum(unif.rayleigh2)))     ct.rayleigh2     <-(chisq.test(obs.rayleigh2,    p = rep(1/length(obs.rayleigh2),    length(obs.rayleigh2))))    $statistic  else  ct.rayleigh2    <-NA
if(length(obs.rayleigh1)     >2 & is.finite(sum(unif.rayleigh1)))     ct.rayleigh1     <-(chisq.test(obs.rayleigh1,    p = rep(1/length(obs.rayleigh1),    length(obs.rayleigh1))))    $statistic  else  ct.rayleigh1    <-NA
if(length(obs.ngamma3)       >2 & is.finite(sum(unif.ngamma3)))       ct.ngamma3       <-(chisq.test(obs.ngamma3,      p = rep(1/length(obs.ngamma3),      length(obs.ngamma3))))      $statistic  else  ct.ngamma3      <-NA
if(length(obs.nweibull3)     >2 & is.finite(sum(unif.nweibull3)))     ct.nweibull3     <-(chisq.test(obs.nweibull3,    p = rep(1/length(obs.nweibull3),    length(obs.nweibull3))))    $statistic  else  ct.nweibull3    <-NA
if(length(obs.nexponential2) >2 & is.finite(sum(unif.nexponential2))) ct.nexponential2 <-(chisq.test(obs.nexponential2,p = rep(1/length(obs.nexponential2),length(obs.nexponential2))))$statistic  else  ct.nexponential2<-NA
if(length(obs.nlognormal3)   >2 & is.finite(sum(unif.nlognormal3)))   ct.nlognormal3   <-(chisq.test(obs.nlognormal3,  p = rep(1/length(obs.nlognormal3),  length(obs.nlognormal3))))  $statistic  else  ct.nlognormal3  <-NA
if(length(obs.nrayleigh2)    >2 & is.finite(sum(unif.nrayleigh2)))    ct.nrayleigh2    <-(chisq.test(obs.nrayleigh2,   p = rep(1/length(obs.nrayleigh2),   length(obs.nrayleigh2))))   $statistic  else  ct.nrayleigh2   <-NA

  CS.Test.Statistic<-c(ct.beta4,ct.beta2,ct.uniform2,ct.uniform1,ct.triangular,ct.gamma3,ct.gamma2,ct.erlang2,ct.weibull3,ct.weibull2,ct.exponential2,ct.exponential1,ct.normal,ct.t3,ct.t1,ct.lognormal3,ct.lognormal2, ct.logistic, ct.rayleigh2, ct.rayleigh1,ct.ngamma3,ct.nweibull3,ct.nexponential2,ct.nlognormal3,ct.nrayleigh2)

# CHI-SQUARE P-VALUE FOR EACH DISTRIBUTION:
if(length(obs.beta4)        >=5 & is.finite(ct.beta4))         P.beta4         <- 1-pchisq(ct.beta4,        (length(obs.beta4)-4-1))         else  P.beta4        <-NA
if(length(obs.beta2)        >=3 & is.finite(ct.beta2))         P.beta2         <- 1-pchisq(ct.beta2,        (length(obs.beta2)-2-1))         else  P.beta2        <-NA
if(length(obs.uniform2)     >=3 & is.finite(ct.uniform2))      P.uniform2      <- 1-pchisq(ct.uniform2,     (length(obs.uniform2)-2-1))      else  P.uniform2     <-NA
if(length(obs.uniform1)     >=2 & is.finite(ct.uniform1))      P.uniform1      <- 1-pchisq(ct.uniform1,     (length(obs.uniform1)-1-1))      else  P.uniform1     <-NA
if(length(obs.triangular)   >=4 & is.finite(ct.triangular))    P.triangular    <- 1-pchisq(ct.triangular,   (length(obs.triangular)-3-1))    else  P.triangular   <-NA
if(length(obs.gamma3)       >=4 & is.finite(ct.gamma3))        P.gamma3        <- 1-pchisq(ct.gamma3,       (length(obs.gamma3)-3-1))        else  P.gamma3       <-NA
if(length(obs.gamma2)       >=3 & is.finite(ct.gamma2))        P.gamma2        <- 1-pchisq(ct.gamma2,       (length(obs.gamma2)-2-1))        else  P.gamma2       <-NA
if(length(obs.erlang2)      >=3 & is.finite(ct.erlang2))       P.erlang2       <- 1-pchisq(ct.erlang2,      (length(obs.erlang2)-2-1))       else  P.erlang2      <-NA
if(length(obs.weibull3)     >=4 & is.finite(ct.weibull3))      P.weibull3      <- 1-pchisq(ct.weibull3,     (length(obs.weibull3)-3-1))      else  P.weibull3     <-NA
if(length(obs.weibull2)     >=3 & is.finite(ct.weibull2))      P.weibull2      <- 1-pchisq(ct.weibull2,     (length(obs.weibull2)-2-1))      else  P.weibull2     <-NA
if(length(obs.exponential2) >=3 & is.finite(ct.exponential2))  P.exponential2  <- 1-pchisq(ct.exponential2, (length(obs.exponential2)-2-1))  else  P.exponential2 <-NA
if(length(obs.exponential1) >=2 & is.finite(ct.exponential1))  P.exponential1  <- 1-pchisq(ct.exponential1, (length(obs.exponential1)-1-1))  else  P.exponential1 <-NA
if(length(obs.normal)       >=3 & is.finite(ct.normal))        P.normal        <- 1-pchisq(ct.normal,       (length(obs.normal)-2-1))        else  P.normal       <-NA
if(length(obs.t3)           >=4 & is.finite(ct.t3))            P.t3            <- 1-pchisq(ct.t3,           (length(obs.t3)-3-1))            else  P.t3           <-NA
if(length(obs.t1)           >=2 & is.finite(ct.t1))            P.t1            <- 1-pchisq(ct.t1,           (length(obs.t1)-1-1))            else  P.t1           <-NA
if(length(obs.lognormal3)   >=4 & is.finite(ct.lognormal3))    P.lognormal3    <- 1-pchisq(ct.lognormal3,   (length(obs.lognormal3)-3-1))    else  P.lognormal3   <-NA
if(length(obs.lognormal2)   >=3 & is.finite(ct.lognormal2))    P.lognormal2    <- 1-pchisq(ct.lognormal2,   (length(obs.lognormal2)-2-1))    else  P.lognormal2   <-NA
if(length(obs.logistic)     >=3 & is.finite(ct.logistic))      P.logistic      <- 1-pchisq(ct.logistic,     (length(obs.logistic)-2-1))      else  P.logistic     <-NA
if(length(obs.rayleigh2)    >=3 & is.finite(ct.rayleigh2))     P.rayleigh2     <- 1-pchisq(ct.rayleigh2,    (length(obs.rayleigh2)-2-1))     else  P.rayleigh2    <-NA
if(length(obs.rayleigh1)    >=2 & is.finite(ct.rayleigh1))     P.rayleigh1     <- 1-pchisq(ct.rayleigh1,    (length(obs.rayleigh1)-1-1))     else  P.rayleigh1    <-NA
if(length(obs.ngamma3)      >=4 & is.finite(ct.ngamma3))       P.ngamma3       <- 1-pchisq(ct.ngamma3,      (length(obs.ngamma3)-3-1))       else  P.ngamma3      <-NA
if(length(obs.nweibull3)    >=4 & is.finite(ct.nweibull3))     P.nweibull3     <- 1-pchisq(ct.nweibull3,    (length(obs.nweibull3)-3-1))     else  P.nweibull3    <-NA
if(length(obs.nexponential2)>=3 & is.finite(ct.nexponential2)) P.nexponential2 <- 1-pchisq(ct.nexponential2,(length(obs.nexponential2)-2-1)) else  P.nexponential2<-NA
if(length(obs.nlognormal3)  >=3 & is.finite(ct.nlognormal3))   P.nlognormal3   <- 1-pchisq(ct.nlognormal3,  (length(obs.nlognormal3)-3-1))   else  P.nlognormal3  <-NA
if(length(obs.nrayleigh2)   >=3 & is.finite(ct.nrayleigh2))    P.nrayleigh2    <- 1-pchisq(ct.nrayleigh2,   (length(obs.nrayleigh2)-2-1))    else  P.nrayleigh2   <-NA

  CS.P.Value<-c(P.beta4,P.beta2,P.uniform2,P.uniform1,P.triangular,P.gamma3,P.gamma2,P.erlang2,P.weibull3,P.weibull2,P.exponential2,P.exponential1,P.normal,P.t3,P.t1,P.lognormal3,P.lognormal2, P.logistic, P.rayleigh2, P.rayleigh1,P.ngamma3,P.nweibull3,P.nexponential2,P.nlognormal3,P.nrayleigh2)
  

# SAMPLE SIZE:
  n <- length(x)  

# CRITICAL VALUES FOR CHI-SQUARE TEST: THE LEVEL OF SIGNIFICANCE (a):
   if   (n<=20000)          a <- 1000/n       #%
   else                     a <- 0.05         #%
   CS.Sig.level<-rep(0.01*a,length(CS.P.Value))

# CHI-SQUARE CRITICAL VALUE FOR EACH DISTRIBUTION:
if(length(obs.beta4)        >=5 & is.finite(ct.beta4))         cv.beta4         <-qchisq((1-0.01*a),(length(obs.beta4)-4-1))         else  cv.beta4        <-NA
if(length(obs.beta2)        >=3 & is.finite(ct.beta2))         cv.beta2         <-qchisq((1-0.01*a),(length(obs.beta2)-2-1))         else  cv.beta2        <-NA
if(length(obs.uniform2)     >=3 & is.finite(ct.uniform2))      cv.uniform2      <-qchisq((1-0.01*a),(length(obs.uniform2)-2-1))      else  cv.uniform2     <-NA
if(length(obs.uniform1)     >=2 & is.finite(ct.uniform1))      cv.uniform1      <-qchisq((1-0.01*a),(length(obs.uniform1)-1-1))      else  cv.uniform1     <-NA
if(length(obs.triangular)   >=4 & is.finite(ct.triangular))    cv.triangular    <-qchisq((1-0.01*a),(length(obs.triangular)-3-1))    else  cv.triangular   <-NA
if(length(obs.gamma3)       >=4 & is.finite(ct.gamma3))        cv.gamma3        <-qchisq((1-0.01*a),(length(obs.gamma3)-3-1))        else  cv.gamma3       <-NA
if(length(obs.gamma2)       >=3 & is.finite(ct.gamma2))        cv.gamma2        <-qchisq((1-0.01*a),(length(obs.gamma2)-2-1))        else  cv.gamma2       <-NA
if(length(obs.erlang2)      >=3 & is.finite(ct.erlang2))       cv.erlang2       <-qchisq((1-0.01*a),(length(obs.erlang2)-2-1))       else  cv.erlang2      <-NA
if(length(obs.weibull3)     >=4 & is.finite(ct.weibull3))      cv.weibull3      <-qchisq((1-0.01*a),(length(obs.weibull3)-3-1))      else  cv.weibull3     <-NA
if(length(obs.weibull2)     >=3 & is.finite(ct.weibull2))      cv.weibull2      <-qchisq((1-0.01*a),(length(obs.weibull2)-2-1))      else  cv.weibull2     <-NA
if(length(obs.exponential2) >=3 & is.finite(ct.exponential2))  cv.exponential2  <-qchisq((1-0.01*a),(length(obs.exponential2)-2-1))  else  cv.exponential2 <-NA
if(length(obs.exponential1) >=2 & is.finite(ct.exponential1))  cv.exponential1  <-qchisq((1-0.01*a),(length(obs.exponential1)-1-1))  else  cv.exponential1 <-NA
if(length(obs.normal)       >=3 & is.finite(ct.normal))        cv.normal        <-qchisq((1-0.01*a),(length(obs.normal)-2-1))        else  cv.normal       <-NA
if(length(obs.t3)           >=4 & is.finite(ct.t3))            cv.t3            <-qchisq((1-0.01*a),(length(obs.t3)-3-1))            else  cv.t3           <-NA
if(length(obs.t1)           >=2 & is.finite(ct.t1))            cv.t1            <-qchisq((1-0.01*a),(length(obs.t1)-1-1))            else  cv.t1           <-NA
if(length(obs.lognormal3)   >=4 & is.finite(ct.lognormal3))    cv.lognormal3    <-qchisq((1-0.01*a),(length(obs.lognormal3)-3-1))    else  cv.lognormal3   <-NA
if(length(obs.lognormal2)   >=3 & is.finite(ct.lognormal2))    cv.lognormal2    <-qchisq((1-0.01*a),(length(obs.lognormal2)-2-1))    else  cv.lognormal2   <-NA
if(length(obs.logistic)     >=3 & is.finite(ct.logistic))      cv.logistic      <-qchisq((1-0.01*a),(length(obs.logistic)-2-1))      else  cv.logistic     <-NA
if(length(obs.rayleigh2)    >=3 & is.finite(ct.rayleigh2))     cv.rayleigh2     <-qchisq((1-0.01*a),(length(obs.rayleigh2)-2-1))     else  cv.rayleigh2    <-NA
if(length(obs.rayleigh1)    >=2 & is.finite(ct.rayleigh1))     cv.rayleigh1     <-qchisq((1-0.01*a),(length(obs.rayleigh1)-1-1))     else  cv.rayleigh1    <-NA
if(length(obs.ngamma3)      >=4 & is.finite(ct.ngamma3))       cv.ngamma3       <-qchisq((1-0.01*a),(length(obs.ngamma3)-3-1))       else  cv.ngamma3      <-NA
if(length(obs.nweibull3)    >=4 & is.finite(ct.nweibull3))     cv.nweibull3     <-qchisq((1-0.01*a),(length(obs.nweibull3)-3-1))     else  cv.nweibull3    <-NA
if(length(obs.nexponential2)>=3 & is.finite(ct.nexponential2)) cv.nexponential2 <-qchisq((1-0.01*a),(length(obs.nexponential2)-2-1)) else  cv.nexponential2<-NA
if(length(obs.nlognormal3)  >=3 & is.finite(ct.nlognormal3))   cv.nlognormal3   <-qchisq((1-0.01*a),(length(obs.nlognormal3)-3-1))   else  cv.nlognormal3  <-NA
if(length(obs.nrayleigh2)   >=3 & is.finite(ct.nrayleigh2))    cv.nrayleigh2    <-qchisq((1-0.01*a),(length(obs.nrayleigh2)-2-1))    else  cv.nrayleigh2   <-NA

  CS.Critical.Value<-c(cv.beta4,cv.beta2,cv.uniform2,cv.uniform1,cv.triangular,cv.gamma3,cv.gamma2,cv.erlang2,cv.weibull3,cv.weibull2,cv.exponential2,cv.exponential1,cv.normal,cv.t3,cv.t1,cv.lognormal3,cv.lognormal2, cv.logistic, cv.rayleigh2, cv.rayleigh1,cv.ngamma3,cv.nweibull3,cv.nexponential2,cv.nlognormal3,cv.nrayleigh2)
  
# SHOW CHI-SQUARE TEST STATISTICS AND CORRESPONDING CRITICAL VALUES IN THE SAME MATRIX
  C.S<-cbind(CS.Test.Statistic,CS.Critical.Value,CS.P.Value,CS.Sig.level)

# CHI-SQUARE TEST DECISION:
   b <- length(CS.Critical.Value)
   CS.Decision <- 0
   for (i in 1:b){
    if    (!is.na(C.S[i,1])){   
          if    (C.S[i,1]<C.S[i,2])     CS.Decision[i] <-"accept"
          else                          CS.Decision[i] <-"reject"
    }
    else
    CS.Decision[i] <-"No fit"
   }
result <-data.frame(C.S,CS.Decision)
rownames(result)<-c("Beta(4)","Beta(2)","Uniform(2)","Uniform(1)","Triangular","Gamma(3)","Gamma(2)","Erlang(2)","Weibull(3)","Weibull(2)","Exponential(2)","Exponential(1)","Normal","Student's-t(3)","Student's-t(1)","Lognormal(3)","Lognormal(2)","Logistic","Rayieigh(2)","Rayleigh(1)","N-Gamma(3)","N-Weibull(3)","N-Exponential(2)","N-Lognormal(3)","N-Rayleigh(2)")
result             
}

cs<-cs.test(x)

# OUR FIT ALL FUNCTION:
fitall<-function(x){

### PARAMETRIC MODEL
#  Find values of AIC and MLEs  
   m <- alldistributions(x)

#  Compute Chi-Square P-values and show whether each model passes or fails Chi-Square Test:
   cs <- cs.test(x)
   CS.P.value  <- cs[,3]
   CS.Decision <- cs[,5]

#  Gather values of AIC, MLEs, Chi-Square Test and names of parameters in the same matrix:
   k <- data.frame(m,CS.P.value,CS.Decision,parnames)

#  Order results in ascending order with respect to AIC values:
   ms <- k[order(k[,1]),]

#  Display all values by omitting "NA"s and display all numbers in 2 digits:
   e <- (ms[!is.na(ms[,1]),])
   a<-data.frame(round(e[,1:5],2),e[,6:10])
   

#  Gather values of AIC, MLEs, Chi-Square Test and Delta AIC in the same matrix.
   f<-a[,1:7]
   #  Calculate delta AIC to see how different the other models are from the the best-fitted model:
   Delta.AIC <- f[,1]-f[1,1]
   ALL <- cbind(f,Delta.AIC)

# If the best fitting distribution passes Chi-Square Test, then recommend it by showing AIC values and parameters. 

     if (ALL[1,7]=="accept") {            
     distribution.name <- rownames(ALL[])[1]
     
     # Display AIC value:
     AIC<-(ALL[1,1])

     #parameters:
     par<-ALL[1,2:5]
     # Display parameters  by omitting "NA"s:
     parameters<-(par[!is.na(par[1,1:4])])
     names<-e[1,8:11]
     par.names<-names[!is.na(names[1,1:4])]
     full.definition.parameters<-paste(par.names,parameters,sep="=",collapse=",")
     
     recommend<-paste("It is recommended to use the", distribution.name,"distribution with the parameters:",full.definition.parameters)
     res<-list("model"=recommend,"AIC"=AIC,"parameters"=parameters,"summary"=ALL)
     res
     }


# Else, recommend resampling with noise by displaying the value of bandwidth.
     else {  
  ### NON-PARAMETRIC MODEL: RESAMPLING with NOISE
    n <- length(x)                   #sample size
    s <- sd(x)                       #standard deviation
    R <- IQR(x)                      #interquartile range
    b <- 1.06*min(s,R/1.34)*n^(-1/5) #bandwidth
    Y <- x[ceiling(runif(10000)*n)]+b*rnorm(10000)

     # Display the value of bandwidth in 3 digits                    
     bandwidth<-round(b,3)
     recommend<-paste("It is recommended to use resampling with noise", "with bandwidth value of", bandwidth)
     res<-list("model"=recommend,"bandwidth"=bandwidth,"summary"=ALL)
     res
     }

return (res)
 }

fitall(x)

